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mass Wilson fermions with non-degenerate fermion masses, as it was used in recent dynamical 
simulations for Nf = 2 + I + I fermion flavors. 
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This article is organized as follows: in the first part we will briefly review the idea of using Multi- 
Step stochastic correction, for more details we refer to our paper [|l|]. We provide a description of 
the Polynomial Hybrid Monte-Carlo algorithm with stochastic correction , which was successfully 
applied in [§] to simulate Nf = 2+ \ + \ quark flavors using the twisted-mass fermion action. In 
the remainder we will summarize some of the results of these simulations, again for more details 
we refer to 

1. Multi-Step stochastic correction 

The idea of Multi-Step stochastic correction is a natural extension of the Two-Step technique 
well established for Multi-Boson algorithms (TSMB-alg.), cf. ^ and references therein. One 
uses multiple correction steps with increasing precision. Our choice to control the precision in each 



step is to use polynomials Pi of increasing order ni (see Sec. 1.4) to approximate the required power 
in the pseudo-fermionic part of the action: 

PiiQ') ^ iQ'r", m^) ^ [{Q^rPi{Q^)---Pi-i{&)\\ a-D 

where Q denotes the Hermitian fermion-matrix (containing either one fermion-flavor or in case 
of the twisted mass formulation one flavor doublet). For example, a = ^2 results in a one flavor 
(twisted-mass: one flavor-doublet) action, a = 1 leads to Nf = 2 (twisted-mass: two degenerate 
flavor-doublets). 

The update-sequence now consists of several nested noisy corrections (cf. [Q| and references 
therein) using the various polynomials.^ As an example we show a sequence using three-step cor- 
rection: 

start ^_^^P2 > ^_^P2 ^ P3 ^ P3 (1.2) 

MstepsPi WistepsPi {NiN2)xPi,N2xP2 

V ' 

^9 steps 

^ V ' 

steps 

Each time a configuration is rejected it will be replaced with the last accepted one by that same 
correction step. 

In the following subsections, we describe applications of this technique to Multi-Boson as well 
as Hybrid Monte-Carlo based algorithms. 

1.1 Multi-Step Multi-Boson algorithm 

The update for the first (Multi-Boson) step using Pi as well as for the gauge-fields are ex- 
actly the same as for the TSMB-algorithm as described in [Q] and applied in we just add 
more correction-steps to the procedure. This changes the numerical cost expressed in units of 
matrix-vector-multiplications (MVMs) given for TSMB in [Q] (for unexplained notation confer that 
publication) to 

cost ' 

= 6{niN^+Nu)+lGFG + Y,{nk + nk)Nck, (1-3) 

k=2 



'For the noisy correction step we also need polynomials P/(x) ~ Pi{x) '/^ of order n; to generate the noise vector 
from a Gaussian distribution. 
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Nck is the number of times the correction step involving P/t is called. In |[T]] we gave an example 
how calling the most expensive correction step less often actually works in a numerical simulation. 

1.2 HMC and PHMC with stochastic correction 

The first update-step does not necessarily have to rely on a Multi-Boson update. It can for in- 
stance be replaced by a standard Hybrid Monte-Carlo( HMC )-update. The advantage of combining 
stochastic corrections with HMC algorithms is that one may use a mass-preconditioned fermion- 
matrix + /^i in the basic update step, and account for this in the stochastic correction steps, 
using there subsequently decreasing /i,- to achieve the "target action" = 0) in the last (N^^) step 

(pLi eR, 1^1 > . . . > Hi > . . . > I^N = 0). 

A severe limitation of the traditional HMC algorithm is, that only an even number of fermionic 
fields can be simulated. To overcome this drawback one has to take a fractional power of the 
squared fermion matrix, which can be done by using an approximation of the function x^", e.g. 
a = 1/2- Commonly, either a polynomial or a rational approximation is used, leading to a so-called 



PHMC or RHMC algorithm, respectively. The PHMC was suggested in m- 12], while a descrip 



tion of the RHMC can be found in [13|. The original proposal of the PHMC algorithm had to face 
the problem, that large orders in the polynomial approximation forbid an efficient application to 
the region of light quark masses. In the remainder of this section we will describe our approach 
to introduce a stochastic correction step to make this algorithm competitive to other algorithms 
currently used. For results on a successful application see Sec. ^. An alternative to the stochas- 
tic correction would be to reweight the generated configurations, this is currently investigated by 



Chiarappa et ai, see [|14|]. 

The HMC algorithm and its variants all start by introducing a fictitious time-scale T and evolve 
the gauge and conjugate momenta fields according to Hamilton's equations of motion. Here we 
will only describe the details concerning the polynomial approximation and correction step in the 
fermionic part of the action. Again, we use a polynomial approximation Pi{Q^) ~ (Q^)^" of order 
ni , where Pi may be evaluated in a recursive scheme (which will be convenient in the Metropolis 
step) or by using the root-representation 

PiiQ^) = c,f[{Q^-n) = cof[{Q-pd fliQ-p*)- (1-4) 

i=\ (=1 i=ni 

The last rearrangement is possible because for wi even (to which we restrict ourselves) the roots 
always appear in complex conjugate pairs. Note the ordering in the two products (first increas- 
ing index, last decreasing index): this assures a better numerical stability (the roots r, itself are 



also ordered to a scheme to minimize the numerical error [|15j|) and allows for a efficient way of 
calculating the fermionic force, as needed in the HMC-step. By introducing the auxiliary fields 

^l''^ = V^^{Q-pi)---iQ-pk), k=l,...,n-l (1.5) 

^ V^oHQ-Pi)---{Q-Pn){Q-p:)---{Q-p*k+2) (1-6) 

and setting ^j*'^ = , the fermionic force can be written as 




(1.7) 
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Here always denotes the pseudo-fermionic fields from the "bosonization" of the fermion matrix 
and D is the derivative with respect to the given indices. Observe, that it is most efficient to first 
calculate all the 0^*^' and store the results. In the actual computation one than has to combine the 

(k) ik) 

actual 02 (which is obtained recursively from the previous one) with one of the prepared 0j - 
fields. Using a better conditioned fermion matrix by applying even-odd-preconditioning is also 
possible, in that case a similar scheme can be written down. 

After a PHMC-step, which starts with generating the pseudo-fermionic fields according to 
= A(G^) ^7 from a Gaussian-distributed vector T] using a polynomial P\{x) ~ Py{x)^^/^ of or- 
der hi and the real conjugate momenta field P^^j (in the adjoint representation) according to the 
contribution dPexp(— /'^/2), a trajectory of length 8t is performed using the Sexton-Weingarten 
integration scheme with multiple time-scales for the gauge and fermionic parts. At the end of 
every trajectory a global accept-reject step is performed, using the polynomial Pi in its recursive 
form to calculate the fermionic part of the energy. The reason for this is to correct for numerical 
errors originating from the finite step size in the integrator. After A^traj such trajectories a stochastic 
correction step is carried out as described in Sec. |l|, using polynomials P2 and P2 in the recursive 
representation. In this way we are able to correct for the (intentionally) poor approximation in the 
PHMC-step, which allows us to prepare trajectories at moderate numerical cost. An approximate 
formula^ for the cost in number of matrix-vector-multiplications (MVMs) can be given as follows: 



cost 

MVM ^ ^^Bini + m) +Mraj 2«b(«i +«i) + ?iB(3 + 2A/^g)(4?ii - 1) 



(1.8) 



where ng is the number of determinant-breakup used (see below) and Nq denotes the number of 
time-steps for the fermionic integration. 

1.3 Determinant-breakup, mass preconditioning 

An important improvement of the Multi-Boson algorithms and the PHMC- as well as the 
RHMC-algorithms is so-called "determinant-breakup" The fractional powers in the approxi- 
mation allow to use ng sets of pseudo-fermion fields with a replaced hy a/ng. This technique was 
proven to improve the simulations for the TSMB- [||-|8]] and the PHMC- [Q] and also recently for 
the RHMC-algorithm 

Other techniques to speed up the simulation rely on preconditioning of the fermion-matrix. 



we already mentioned mass-preconditioning briefly in Sec. |1.2| . Details on how to use this in 
conjunction with multiple correction steps are given in 

1.4 Polynomial approximation 

There are several different methods available for the required approximations in the MSMB 
and HMC-variant algorithms. Most commonly used are polynomial or rational approximations 
of some order n. Usually for a given precision the rational approximation requires a lower order 



[13 1, but here we want to stress the fact that while the order of a polynomial determines the cost in 



matrix-vector-multiplications directly, for the rational function one has to invert the fermion matrix 



Here we neglected the cost from calculating the force from the gauge and conjugate moment parts and made the 



assumption that the tensor-like multiplication in Eq. (1.7) counts as one MVM. 



4 



Multi-Step stochastic correction in dynamical fermion updating algorithms 



Enno E. Scholz 



0.02 
0.015 

0.01 
0.005 


-0.005 
-0.01 

-0.015 
-0.02 



rel. deviation of 
n■^ = 60 
a. = 0.5 



E 
1 

'ii: 


= 5e-5 
= 4.0 

i ft A n ft A A 




ro 

10.2 
-0.4 
-0.6 
-0.8 




1 



0.01 
0.005 


-0.005 
-0.01 




rei. deviation of P2P3 
ng = 800 
a = 0.5 
E = 5e-5 




Figure 1: Example for polynomials with increasing order to approximate x -1/2 in the interval [5 • 10-^4.0] 
as used in a MSMB-update [|l]], shown is the relative deviation 



at least once, which is the main contribution to the cost and is itself a polynomial approximation, 
too. Therefore we prefer to take a direct polynomial ansatz.-' 

Another question arises, namely in which way the error of the approximation is minimized. In 
principle, defining any norm and calculate the deviation according to it is a legitimate procedure. 
Commonly used are either the L2- or the Lt„-norm. The first minimizes the quadratic deviation, 
while the latter minimizes the maximal deviation. It has been shown that the first one is better 



suited for dynamical fermion simulations []18|], since it leads to smaller deviations in the bulk of the 
approximation-interval at the cost of a higher deviation at the boundaries. The Loo-norm leads to an 
uniformly distributed deviation, which is disadvantageous because most of the eigenvalues of the 
fermion matrix are in the bulk and only a few ones are close to the lower boundary. Examples of 
the L2-optimized polynomial approximations as used in a three-step multi-boson correction update 
[|l|] are shown in Fig. |l|. For the generation of the coefficients we refer to [15, |l^, 20]. 



2. Twisted mass-simulations with Nf = 2 + 1 + 1 

As part of the efforts of the European Twisted Mass-Collaboration (ETM-ColL), the PHMC- 
algorithm with one stochastic correction step as introduced above has been applied to simulate two 
doublets of twisted mass quarks [Q]. In the twisted mass action [21] an additional mass term is 
added, which results in a lower bound for the eigenvalue spectrum and, if the standard mass term 
is properly tuned, in an automatically ^(a)-improved fermion action. The first doublet contains 
two mass-degenerate Ught quark flavors (up and down), whereas in the second doublet a mass- 



splitting term as introduced by Frezzotti and Rossi [ ]22| ] has been added. Therefore, the quarks of 
the second doublet can be identified as the charm and strange flavors, having different masses. In 
the mass-split case we have to use an algorithm which is capable of simulating a single fermion 
flavor. In principle it would have been possible to use the HMC as in [ p3| ] for the first doublet and 
the PHMC only for the second one. For simplicity we have chosen to use the PHMC-algorithm for 
both doublets and the results showed that the PHMC with stochastic correction works fine for light 
quarks, too. 



^Note, that although the polynomial degrees of the first approximation steps are fixed, the degree of the last (and 
most expensive) step can be adaptive because of the recursive evaluation. 
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Figure 2: Squared pion- and kaon-masses and their ratio from 12^ x 24 and 16^ x 32 lattices, for details see 
text and [0] 



We performed simulations at two lattice spacings on a fixed physical volume of aL ~ 2.4fm 
(a « 0.20fm at x T = 12^ x 24 and a w 0.15fm at 16^ x 32) using the Symanzik tree-level 



improved gauge action [|24|]. In both cases the twisted mass parameters where kept fixed in physical 
units and the untwisted mass parameter was varied. On the coarse lattice the numerical cost"* to 
generate a new configuration varied between 28.2 and 46.5 (in thousands of MVMs) for the heaviest 
and lightest pion mass. On the finer lattice we observed similar numbers between 17.9 and 44.2. 
For more details see [^. That allowed us to study the phase structure and determine the critical 
value for the untwisted mass parameter, which is important to obtain the ^ (a) -improvement. 

The volume scaling of the algorithm is very good: keeping the action parameters constant, 
for instance, between 12^ x 24 and 24-^ x 48 lattices, one can keep the polynomial degrees n2, «2 
constant and increase {n\+h\) only by about 10% for the same lower bound of the approximation 
intervall. Even this moderate increase can be compensated by increasing the lower bound which is 
possible due to the diminished fluctuations of the smallest eigenvalues. 

2.1 Results 

Here we summarize the results, details can be found in The most important finding con- 
cerns the phase structure: at the coarser lattice spacing a strong metastability occurs (as was already 



found for Nf = 2 Wilson fermions, cf. ||8[ g5|, |2^ and references therein), leading to a minimal 



pion mass of approx. 670MeV. At the finer lattice spacing the metastability could not be ob- 
served. Actually, since the lattice volume was fixed, we were not able to distinguish there between 
a phase-transition or a cross-over scenario. Anyway, without a metastability a lightest pion mass 
of 450MeV has been achieved. This implies that on a 24^ x 48 lattice with a ~ O.lfm a pion mass 
of about 300MeV can probably be simulated. Fig. § shows the squared pion- and kaon-masses and 
the ratio (m^j/mTf)^ as functions of the untwisted quark mass. The latter plot also contains chiral 
perturbation theory guided fits to the data. Remarkably, the minimum value of these fits is close to 
the physical point {niji/niKf' — 0.08. 

To conclude, dynamical u, d, c, and s quarks can be simulated in the Frezzotti-Rossi twisted 
mass formulation with a moderate tuning effort. The PHMC algorithm is working fine for both 
light (m, d) and heavier (c, s) quarks. 



^To compare these numbers with those from other groups: IMVM ~ 2688 • flflop, where Q. = xT is the lattice- 
volume. 
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